# set working directory
setwd('~/replication_files')

# function to load and install packages
load_packages<-function(packages){
  sapply(packages,function(x){
    f<-require(x,character.only=T)
    if(!f) install.packages(x); f<-require(x,character.only=T)
    return(f)},USE.NAMES=T)
}

# load packages
load_packages(c('data.table','ggplot2','gtools','plotrix','gridExtra','diagis','stringr'))

# function to extract legend from ggplot object
g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

FinDT<-readRDS('Figure1.rds')
alloc.col(FinDT,ncol(FinDT)+20)

# create weighted means of pro and antivax page views using survey weights by tercile of vaccine attitudes
# in addition, create the lower and upper 95% confidence values
# first do this for the binary measure (transformed to percentage), then the count measure
FinDT[,provax_bin_terc_pct:=100*weighted.mean(provax_bin,weight),by='vax_tercile']
FinDT[,provax_bin_terc_pct_hi:=100*(weighted.mean(provax_bin,weight)+qnorm(.975)*weighted_se(provax_bin,weight)),by='vax_tercile']
FinDT[,provax_bin_terc_pct_lo:=100*(weighted.mean(provax_bin,weight)-qnorm(.975)*weighted_se(provax_bin,weight)),by='vax_tercile']
FinDT[,antivax_bin_terc_pct:=100*weighted.mean(antivax_bin,weight),by='vax_tercile']
FinDT[,antivax_bin_terc_pct_hi:=100*(weighted.mean(antivax_bin,weight)+qnorm(.975)*weighted_se(antivax_bin,weight)),by='vax_tercile']
FinDT[,antivax_bin_terc_pct_lo:=100*(weighted.mean(antivax_bin,weight)-qnorm(.975)*weighted_se(antivax_bin,weight)),by='vax_tercile']

FinDT[,provax_count_terc:=weighted.mean(provax_count,weight),by='vax_tercile']
FinDT[,provax_count_terc_hi:=weighted.mean(provax_count,weight)+qnorm(.975)*weighted_se(provax_count,weight),by='vax_tercile']
FinDT[,provax_count_terc_lo:=weighted.mean(provax_count,weight)-qnorm(.975)*weighted_se(provax_count,weight),by='vax_tercile']
FinDT[,antivax_count_terc:=weighted.mean(antivax_count,weight),by='vax_tercile']
FinDT[,antivax_count_terc_hi:=weighted.mean(antivax_count,weight)+qnorm(.975)*weighted_se(antivax_count,weight),by='vax_tercile']
FinDT[,antivax_count_terc_lo:=weighted.mean(antivax_count,weight)-qnorm(.975)*weighted_se(antivax_count,weight),by='vax_tercile']

# get the unique values of these variables
(FinDT<-unique(FinDT[,list(vax_tercile,provax_count_terc,antivax_count_terc,provax_bin_terc_pct,antivax_bin_terc_pct,
                           provax_count_terc_hi,provax_count_terc_lo,antivax_count_terc_hi,antivax_count_terc_lo,
                           provax_bin_terc_pct_hi,provax_bin_terc_pct_lo,antivax_bin_terc_pct_hi,antivax_bin_terc_pct_lo)]))

# melt the data.table to create one with the point estimates and two more with the lower and upper bounds
(binDT<-melt(FinDT,id.vars='vax_tercile',measure.vars=c('provax_bin_terc_pct','antivax_bin_terc_pct')))
(binDTlo<-melt(FinDT,id.vars='vax_tercile',measure.vars=c('provax_bin_terc_pct_lo','antivax_bin_terc_pct_lo'),value.name='lo'))
(binDThi<-melt(FinDT,id.vars='vax_tercile',measure.vars=c('provax_bin_terc_pct_hi','antivax_bin_terc_pct_hi'),value.name='hi'))

# add the lower and upper bound values to the data.table
set(binDT,NULL,'lo',binDTlo[['lo']])
set(binDT,NULL,'hi',binDThi[['hi']])
# create factor of vaccine attitude tercile, add labels
set(binDT,NULL,'vax_tercile_fac',factor(binDT[['vax_tercile']],levels=1:3,labels=c('Least favorable','Somewhat favorable','Most favorable')))
# create factor of the type of page, add label
set(binDT,NULL,'vax_type_fac',factor(rep(1:2,each=3),levels=1:2,labels=c('Not vaccine-skeptical','Vaccine-skeptical')))

# do what we did above for the binary data (the percentages) for the count data:
(countDT<-melt(FinDT,id.vars='vax_tercile',measure.vars=c('provax_count_terc','antivax_count_terc')))
(countDTlo<-melt(FinDT,id.vars='vax_tercile',measure.vars=c('provax_count_terc_lo','antivax_count_terc_lo'),value.name='lo'))
(countDThi<-melt(FinDT,id.vars='vax_tercile',measure.vars=c('provax_count_terc_hi','antivax_count_terc_hi'),value.name='hi'))

set(countDT,NULL,'lo',countDTlo[['lo']])
set(countDT,NULL,'hi',countDThi[['hi']])
set(countDT,NULL,'vax_tercile_fac',factor(countDT[['vax_tercile']],levels=1:3,labels=c('Least favorable','Somewhat favorable','Most favorable')))
set(countDT,NULL,'vax_type_fac',factor(rep(1:2,each=3),levels=1:2,labels=c('Not vaccine-skeptical','Vaccine-skeptical')))

# create a plot to get the legend
p1.leg<-
  ggplot(data=binDT,mapping=aes(x=vax_tercile_fac,y=value,group=vax_type_fac))+
  geom_col(mapping=aes(color=vax_type_fac,fill=vax_type_fac),position='dodge')+
  geom_linerange(aes(ymin=lo, ymax=hi),position=position_dodge(width=.9),color='gray80',size=.6)+
  scale_color_manual(values=c('gray60','black'))+
  scale_fill_manual(values=c('gray60','black'))+
  guides(color=guide_legend(title='Webpage type'),fill=guide_legend(title='Webpage type'))+
  theme_bw()+theme(axis.ticks=element_blank(),panel.border=element_blank(),panel.grid.major.x=element_blank(),
                   panel.grid.minor=element_blank(),legend.position='bottom',axis.text=element_text(color='black'))

# create the plot for the percentage data
p1<-
  ggplot(data=binDT,mapping=aes(x=vax_tercile_fac,y=value,group=vax_type_fac))+
  geom_col(mapping=aes(color=vax_type_fac,fill=vax_type_fac),position='dodge')+
  geom_linerange(aes(ymin=lo, ymax=hi),position=position_dodge(width=.9),color='gray80',size=.6)+
  scale_color_manual(values=c('gray60','black'))+
  scale_fill_manual(values=c('gray60','black'))+
  guides(color=F,fill=F)+
  labs(x='Vaccine attitude',y=NULL,title='Visited at least one page')+
  scale_y_continuous(labels = function(x) paste0(x,"%"))+
  theme_bw()+theme(axis.ticks=element_blank(),panel.border=element_blank(),panel.grid.major.x=element_blank(),plot.title=element_text(hjust=0.5),
                   panel.grid.minor=element_blank(),axis.text=element_text(color='black'),axis.text.x=element_text(margin=margin(-10,0,5,0)))

# create a plot for the count data
p2<-
  ggplot(data=countDT,mapping=aes(x=vax_tercile_fac,y=value,group=vax_type_fac))+
  geom_col(mapping=aes(color=vax_type_fac,fill=vax_type_fac),position='dodge')+
  geom_linerange(mapping=aes(ymin=lo, ymax=hi),position=position_dodge(width=.9),color='gray80',size=.6)+
  scale_color_manual(values=c('gray60','black'))+
  scale_fill_manual(values=c('gray60','black'))+
  guides(color=F,fill=F)+
  coord_cartesian(ylim=c(0,1.18))+
  scale_y_continuous(breaks=seq(0,1.2,.3))+
  labs(x='Vaccine attitude',y=NULL,title='Mean number of pages visited')+
  theme_bw()+theme(axis.ticks=element_blank(),panel.border=element_blank(),panel.grid.major.x=element_blank(),plot.title=element_text(hjust=0.5),
                   panel.grid.minor=element_blank(),axis.text=element_text(color='black'),axis.text.x=element_text(margin=margin(-10,0,5,0)))

# arrange plots and legend together
p<-arrangeGrob(arrangeGrob(p1,p2,nrow=1),g_legend(p1.leg),ncol=1,heights=c(1,.08))

# save plot
ggsave('Figure1.pdf',p,width=8,height=6,units='in',dpi=300)
